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Abstract 

We analyze the dynamical evolution of binary stars that interact with a static background of single stars in 
the environment of a massive black hole (MBH). All stars are considered to be single mass, Newtonian point 
particles. We follow the evolution of the energy E and angular momentum J of the center of mass of the 
binaries with respect to the MBH, as well as their internal semi-major axis a, using a Monte Carlo method. 
For a system like the Galactic center, the main conclusions are the following: (1) The binary fraction can be of 
the order of a few percent outside 0.1 pc, but decreases quickly closer to the MBH. (2) Within ~0.1 pc, binaries 
can only exist on eccentric orbits with apocenters much further away from the MBH. (3) Far away from the 
MBH, loss-cone effects are the dominant mechanism that disrupts binaries with internal velocities close to 
the velocity dispersion. Closer to the MBH, three-body encounters are more effective in disrupting binaries. 
(4) The rate at which hard binaries become tighter is usually less than the rate at which a binary diffuses to 
orbits that are more bound to the MBH. (5) Binaries are typically disrupted before they experience an exchange 
interaction; as a result, the number of exchanges is less than one would estimate from a simple ''nva estimate". 
We give applications of our results to the formation of X-ray binaries near MBHs and to the production rates 
of hyper-velocity stars by intermediate mass MBHs. 

Subject headings: black hole physics — stellar dynamics — Galaxy: center — binaries: general 



1. EVITRODUCTION 

The presence of a massive black hole (MBH) in a stellar 
system has several interesting implications for stellar dynam- 
ics in those systems (see e.g. Alexander 2003). Black holes 
provide the deepest gravitational potential wells in nature. As 
a result, MBHs give rise to a strong tidal field, and extremely 
high velocity dispersions which diverge as Udisp « \f^Jr. 
At the same time, the stellar density in these systems also di- 
verges towards the center In spite of the high velocities, the 
rate of encounters between stars can therefore be sufficiently 
high for the system to exhibit relaxational evolution within the 
age of the Universe. In this paper we explore the dynamical 
evolution of binary stars under the conditions of a very hot, 
dense system with a tidal sink. 

1.1. Collisional stellar dynamics near a massive black hole 

The relaxation time is for a single mass system convention- 
ally defined as 



n{r) cx r-'^l^- N{E) cx E-^'^, 



(2) 



where r is the radial distance from the MBH, E — GM,/r — 
v'^ /2 the negative specific energy with respect to the MBH for 
a star with velocity v (hereafter "energy"); and N{E)dE is the 
number of stars in the interval' {E, E + dE). Extensive theo- 
retical work has confirmed the existence of a Bahcall & Woij 
d 19761) pr ofile under a wide variety of assumptions and meth- 
ods, (e.g. IShai)ir o & Marchan li[T978l:ICohn & KulsrudlfTOTSl: 
Frejtag & Bend 2002; Baum gardt et alj l2004al; iPreto et all 
|2004|) . including iV-body simulations with primordial bi- 
naries (Trentietal. 2007). There is now also observa- 
tional support of a diverg ing density profile th at is approxi- 



1976 ) in the Galac- 
200aiSchoderetan 



(Gm)2n In A' 



(1) 



where m is the mass of a star, n is the number density, 
'^disp is the velocity dispersion, and In A is the Coulomb 
logarithm ( S pitzer & Hart 1 9711 iBinnev & Tremai ne 200^)- 
This time scale is a measure for the time it takes for a sys- 
tem to evolve as a result of the sum of many weak, un- 
corrected two-body interactions. The resulting density pro- 
file of single ma ss star s close to a MBH was first found by 
lahcall & Wolf ( fT976l) . They solved the one-dimensional 
(only energy) Fokker-Planck equations under the approxima- 
tion that the potential of the system is given by the MBH. In 
steady state, the stellar profile is given by 
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mately co nsistent with | Bahcall & Wolf! 
tic ce nter (lAlexandeij|1999HGenzel et alj 
l2007h . It is important to note, however, that the interpre- 
tation of relaxational stellar dynamics in the Galactic cen- 
ter is complicated by the presence of many luminous young 
stars, which domina te the light, but are certainly not re- 
laxed (e.g. [L evin & Beloborodovll2003l: iPaumard et alll2006l: 
iBartko et al.ll200 8). and giants, which may be affected by hy- 
drodynamical collisions (Alexander 1999). For a comprehen- 
sive review of stellar phenomena and observations thereof in 
the Galactic center, see i Alexander! (i2005h . 

From equations ([H |2|l, it can be seen that for a 
iBahcall & Wolj (11976') profile, the relaxation time de- 
creases as r^/'^ as r decreases, in spite of the diverging 
velocity dispersion. Many systems, perhaps most notably 
galaxies as a whole, have relaxation times that are many 
orders of magnitude larger than the age of the Universe, 
and the effect of the diffusion of stars as a result of two- 
body interactions is entirely negligible. However, close 
to MBHs of masses M, < 10'' M©, the relaxation time 
is smaller than the age of the Universe. Such systems 

' We will use the symbol N for several functions, with the argument im- 
plying its meaning. 
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are known as collisional, and give rise to a variety of 
phenomena that are absent for systems with much larger 
relaxation times. These phenomena include the formation 
of the B ahc all & Wolf (1976) cusp; m ass se gre gation (e.g. 
Bahcall & Wolf 1977; Murphv et al. 1991; Freitag et aij 
2006; Hpcman & AlexMder . 2006b; O'Leary et al. "2008t 



Alexander & Hopm anllioost iKeshet etal J 120091); loss-cone 
dynam i cs (e.g. |Frank&Rees 119761; i Lightman & Shapird 
T977t iB^all&Wolfl 119771; iCohn & Kulsrud 1978; 



Sver&Ul meii 119991; iMagorrianl Tremainel 119991: 



Amaro-Seoane et alj l2004l for more details see below); 



inspi ral of stars due to gravitational waves or t i dal dissipation 
(e.g. 'Hils & Be ndal[T99llFreit ag"2001'. "2003^, 'Hopman et al.' 
|2()04; Hopman '&Alexanded 12005; Baumgardt et al. 2005); 
and evolution of binaries, which is the subject of this paper. 

In summary, the dynamics in the close vicinity of a MBH 
is characterized by very high velocity dispersions; very high 
densities; and short relaxation times. These three aspects lead 
to an extreme environment for stellar processes. 

1.2. Binaries 

Binary stars play a major role in many branches of astron- 
omy. Interest in binaries from a dynamical perspective has 
greatl y increased startin g with the seminal papers by Heggie 
(119751) and Hillsj (Il975h . From a dynamical perspective bina- 
ries have a dual role. 

When (single) stars interact with binaries, the internal prop- 
erties of the binaries change. The binary can become tighter 
or less tight, it can dissolve, and it can exchange one of its 
members with an incoming star It is therefore interesting to 
study the evolution of the binaries themselves as they interact 
with their environment. 

In some cases, binaries can also have a dominating role on 
the environment around them. In the point particle approxi- 
mation, a binary can be infinitely tight, and hence provide an 
infinite source of energy. As a result, binari es appear to play 
a major role in stopping core-cqllapse (e.g. lGoodman & Hua 
fT989[lGao et alJll991tlHutetanfT992h . 

In the environment of a MBH that is of interest here, bina- 
ries do not dominate the energetics, since stars can become 
very bound to the MBH, and hence the energies of the sin- 
gle stars also become naturally very large. This conclusion is 
confirmed by iV-body simulations, which show that binaries 
do not affect the global evolution of a stellar system with a 
MBH dTrenti et"al.l 120071) . It is interesting, however, to con- 
sider how they react to their environment. Interactions with 
single stars and with the MBH can modify their distribution 
and composition. 

When all stars are of single mass m, which is what is as- 
sumed throughout this paper, there is a natural dimensionless 
number 



^disp 



Gm 
2^ 



(3) 



that characterizes a binary with semi-major axis a moving 
through a sea of stars that is characterized by a velocity disper- 
sion Udisp- We use the symbol e = Gm? / (2a) for the negative 
internal energy of the binary (hereafter "energy"; note that in 
contrast to the energy E with respect to the MBH, this energy 
is not specific). 

Binaries with ^ ^ 1 are called "soft". It is established 
by extensive theoretical and numerical work that on aver- 
age, such binaries will become even softer, so that ^ 0, 



until they eventually dissolve (e.g. lHeggie|[T97l lHutlll983l; 
lHut&Bahc"aiilll983h . Binaries with ^ » 1 are called 
"hard". In interactions with their environment, such bina- 
ries tend to become even harder (e.g. Heggie 1975; Hut 1993|; 
iHeggie et al.l[T996 ). 

From a dynamical perspective, the evolution of binaries 
near a MBH is interesting for several reasons. 

As outlined above, the close environment near a MBH is 
unique in that it has a diverging density profile, a diverging 
velocity field, and a sufficiently short relaxation time to allow 
for dynamical evolution. Typically, the rate at which binaries 
interact with single stars is of the order of the relaxation rate 
divided by the Coulomb logarithm (see ^for a more quanti- 
tative comparison). A short relaxation time therefore implies 
a high single-binary interaction rate. The fact that the velocity 
dispersion diverges as the radius from the MBH r decreases, 
implies that a binary that is hard at large r, may well be soft 
at smaller r. As this binary diffuses inwards due to two-body 
interactions, it therefore first hardens due to interactions with 
the cold stars around it. But if it does not do so at a sufficiently 
high rate, it will at smaller r be surrounded by very hot stars, 
soften and ionize at smaller r. 

Another aspect that is unique for the system considered in 
this paper, is that at small r the tidal field from the MBH can 
be much larger than the forces that bind the binary compo- 
nents; tidal disruptions therefore form another channel of bi- 
nary destruction (©. 

While the distribution of binary stars is an interesting the- 
oretical problem in stellar dynamics, it also has a number of 
direct observational implications. We will not investigate all 
of these phenomena in this paper, but list them as a motivation 
of the interest in binary dynamics near MBHs. 

X-ray binaries — Within the inner pc of the Galactic cen- 
ter, there is a n over- abundance of X-ray binaries compared 
to the bulge (iMuno et al.l |2005|) . It has long been known 
that globular clusters, which also have small relaxation times, 
contain X-ray binaries in numbers per unit mass that are or- 
ders of magnitude l arger than the rest of the bulge (see e.g. 
iHeggie & Hutll20()3l for a review), and it is thought that col- 
lisional dynamics plays a decisive role in this overabundance 
(Pooley et al. 2003). In this paper the number of stars that 
has experienced an exchange interaction is followed, and from 
this the number of compact remnants present in binaries is es- 
timated. The results are discussed in ^ 

Ultra-luminous X-ray sources — When massive binaries are 
tidally disrupted by an intermediate mass black hole (IMBH) 
of M, = 10"^ — 10^ A/0, the captured star may after ex- 
pansion feed the black hole. This possibility was proposed 
as a possible expla nation for ultra-luminous X-ray sources 
dBlecha et al.ll2006l) . 

Pulsars — Tidal disruption of binaries with compact rem- 
nants may lead to the presence of neutron stars, and perhaps 
millisecond radio pulsars very close to IMBHs. Detection of 
such a pulsar so close to an IMBH w ould give strong ev- 
idence for the presence of the IMBH (iPfahl & Loebl 120041; 
Pfahl 2005; Patruno et al. 2005) 

Hyper velocity stars — If a binary is disrupted in the tidal 
field of a MBH, one of the components can be ejected at a 
very large speed (^ 1 0'^ kms^^) from the Galaxy (Hills 198^ 
lYu & T remaine'2003'). These hyper-velocity stars are now ob- 
served (Brown et al. 2005, 2006, 2007; Brown 2008). In |9] 
we comment on the possibility of ejection of hyper-velocity 
stars from an intermediate mass black hole in the Large Mag- 
ellanic Cloud. 
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S-stars — The energy carried by a hyper-velocity star 
escaping from the energy is drawn from the other com- 
ponent of the original binary, which ends up on a close 
orbit to the MBH. It has been proposed that the clus- 
ter of young, randomly oriented stars in the Galactic 
center, known as the S-stars originated in this process 
fcould & Guillen 2003; Perets et al. 20071 iLockmann etall 
l2008c,Madigan et al.,i2008^.Peretsi,2009i) . 

Gravitational waves: LISA — If one of the components 
of a binary is a compact remnant, and it is on a tight or- 
bit around the MBH after the binary is tidally disrupted, it 
may later spiral in to the MBH while emitting gravitational 
waves. This process may be of considerable importance for 
gravitational wave detection of extreme mass ratio inspiral 
sources by the Laser Interferometer Space Antenna (LISA) 
jMilleret al. 2005). 

Gravitational wav es: LIGO — It was shown by 
lO'Learv et al.l (l2008b that due to the extreme stellar densi- 
ties, occasionally stellar black hole binaries form in the close 
vicinity of the MBH. If such binaries avoid subsequent tidal 
disruption and ionization, and they spiral in within the age of 
the age of the Universe, they may be an important source of 
short wave length gravitational waves that can be observed by 
advanced LIGO. 

Blue stragglers — Mergers of the two components in a 
binary star may lead to rejuvenation of a binary. Such rejuve- 
nat ed stars are known as " blue stragglers" in globular clusters; 
see iHeggie & Hul (l2003h for a review. This phenomena may 
also be of r elevance in th e formation of some of the hyperve- 
locity stars (lPeretsll2008l) . 

1.3. Method 

Several different techniques have been used to study stel- 
lar dynamics in the vicinity of MBHs. The conceptually most 
straightforward method is to integrate the system by di rect A^- 
body simulations (e.g., Baumgardt et al. 2004a b; Pret o et al] 
2004]). A s t udy in cluding primordial binaries was made by 
Trenti et aP (l2007l) . The main disadvantage of direct iV-body 
simulations is that they are limited to an unreahstically small 
number of particles; the largest study bv iTrenti et alj (l2007h 
had 16384 particles, whereas a minimum of several millions 
of particles is required to model the GC. In a situation where 
collisional dynamics is driving the evolution, the outcome 
may depend crucially on the particle number, so that straight- 
forward scaling of the results is not always possible. In ad- 
dition, the finite radius of the stars provides a new length 
scale; as a result, the binary distribution and disruption rate 
near M, = IO^A/q IMBHs and M, = 3 x W^Mq MBHs is 
quite different. 

Other methods that have been used to study stellar dynam- 
ics near MBHs rely explicitly or implicitly on the Fokker- 
Planck approximation. Statistical approaches like solving 
Boltzmann equations or Fokker-Planck equations (mostly us- 
ing a Monte Carlo [MC] approach) have also been used fre- 
quently in the past to find the evolution of the semi-major 
axis of binaries (e.g . Gao et al. 1991; Goodman & Hut 1993; 
llvanovaet a l7 2005; Blecha et al.ll2006l) 

ExpUcit integration of the Fokker-Planck equations in 
one dimension (just the energies of the stars with re- 
spect to the M BH) were made first in the pioneering 
work of iMlcaLl & Wolf (1976, 1977), and lat er by sev- 
eral o ther authors (e.g. [M urphv et al. 1991; Preto etal. 
2004t iHopman & AlexandeJT2006a.b; O'Leary et al.. ,2008; 
Alexander & Hopmanll2008l: iKeshet et"al.ll2009[) . There has 



also been a two dimensional study bv lCohn & Kulsrudl(ll978h 
that followed both the energies and angular momenta of the 
stars with respect to the MBH using direct integration. Di- 
rect integration of the Fokker-Planck equations, along with 
gaseous models ( A maro-Seoane et a l. 2004) is probably the 
fastest method to study the collisional effects of interest here. 
The disadvantage of these methods is that it is hard to imple- 
ment numerically in an efficient way, and that it is harder to 
adapt to include additional physics than MC methods. 

Fokker-Planck equations can also be solved by a MC 
method. In this method, the distribution function (DF) is rep- 
resented by a realization with an arbitrary number of parti- 
cles, that evolves by a combination of deterministic (drift) 
and random steps (diffusion). In recent years, extensive use 
of the lHeno n (1973) implementation of this method has been 
made (e.g. .Freitag & Benz 2001, 2002; Giirkan et al. 200j; 
iFreitag et al.l 120061) . These papers describe this method in 
much detail, and also compare it to other methods. Here, 
we use a different, simpler mod el here that is bas ed on 
'Shapi ro & Marchanj fl978 ) (see Fr eitag & Be"n3l200lL for a 
discussion of the differences between the methods). This 
model is easier to implement, it is suitable for the problem 
at hand, and it keeps a structure that is closer to the origi- 
nal Fokker-Planck equation. Due to easily implemented in- 
dividual time steps, it also allows an accurate treatment of 
stars on very ecc entric orbits, that coul d become tidally dis- 
rupted. The Sha piro & Marchanj (Il978l) method, and the cur- 
rent implementation, is explained in some detail in ^ We 
extend the two-dimensional (energy and angular momentum) 
IShapiro & Marchanj (Il978l) model to three dimensions, the 
third dimension being the internal semi-major axis of the bi- 
nary. 

1.4. Plan 

This paper is organized as follows. In ^ several results 
about single mass stellar systems near MBHs are recapitu- 
lated, and systems with different MBH masses are related 
through the relation between MBH mass and velocity disper- 
sion of the host galaxy. Some of the more important assump- 
tions used in this paper are summarized. In ^we discuss how 
we will treat the evolution of the center of mass of the binaries 
with respect to the MBH as these diffuse through energy and 
angular momentum space. In ^the evolution of the internal 
semi-major axis of the binaries is described; the case of de- 
struction of binaries is treated separately in ^ The dynamics 
in §ffi]-|5]give rise to several time-scales; in ^these time- 
scales are compared, leading to a qualitative overview of the 
dynamics, including a parametrized analytical model for the 
binary disruption rates. In ^we describe the MC scheme, in 
which the dynamics is implemented. In ^ a number of dif- 
ferent cases are studied in detail and the results are presented. 
In ^ the paper is summarized and its results are discussed. 
We apply our results to the formation rate of X-ray binaries 
in the Galactic center, and to the possibility of the ejection of 
hyper-velocity stars by the Large Magellanic Cloud. 

2. SCALING RELATIONS OF A SMGLE MASS CUSP 

Throughout this paper, it is assumed that binary stars inter- 
act only with a static background of single stars, and with a 
MBH. All stars are of Solar type. An important simplification 
that directly follows from this is that binary-binary interac- 
tions are neglected, and as a result that the evolution is linear, 
i.e., the diffusion coefficients in the Fokker-Planck equation 
are independent of the DF of the binaries. Because of the high 
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velocity dispersion near MBHs, only very tight binaries can 
survive. The binary fraction near MBHs, even for marginally 
bound orbits, is therefore expected to be of order of only a few 
percent. This justifies neglecting mutual binary interactions to 
first order 

The distribution of single stars of equal mass m a r ound a 
MBH of mass M. was first found bv lBahcall &"Woig (fT976h . 
who derived and solved the Fokker-Planck equations describ- 
ing their evolution in energy space. The shape of the resulting 
distribution within the radius of influence of the MBH is in- 
dependent of any of the parameters (although it is of course 
subject to a number of assumptions, such as that M, ^ m). 
Here we recapitulate some of their results. To scale systems 
with different MBH masses, the relation between the mass of 
MBHs, A/,, an d the veloc ity dispersion vq of their host bulge 
or galaxy (.Ferrarese & Me rritt .2000; .Gebhardt et al.. .20001: 
iTremaine et al.ll2002l) is used. 



M, = 1.3 X lO^Mg 



vo 



200 km s- 



(4) 



Both the exponent and the pre-factor have small uncertainties, 
which are ignored here. Although the relation (HI has only 
been observed for masses Af, > IO^A/q, we assume that it 
can be extrapolated to lower mass MBHs and IMBHs. 

The relation (jUl can then be exploited to scale several quan- 
tities with MBH mass. The velocity dispersion for unbound 
orbits is 



Vo = 78 kms 



Af. 



1/4 



^3 X 10^ Mq^ 
Unbound stars that come closer than a distance 

1/2 



GM, 

To = — — = 2.1 pc 



A/. 



3 X lO^AU 



(5) 



(6) 



are strongly influenced by the MBH; tq is known as the ra- 
dius of influence. Our analysis will be restricted to distances 
less than tq. The amount of mass in stars within the radius 
of influence is usually comparable to the MBH mass. This 
implies that the density number density no can be scaled with 
the observed density in the Galactic center. 



no = 4.5 X lO'^pc" 



M. 



-1/2 



(7) 



^3 X 106 AfQ^ 

jGenzel et al.l 12003b . where it was assumed that stars are of 
Solar mass. 
The period at the radius of influence is 

.3 \l/2 



" \GM, 

= 1.7 X lO^r 



M. 



3 X 106 Af^ 



1/4 



(8) 



A time-scale of paramount importance for the purpose of this 
paper is the two-body relaxation time (see Eq. [TJ. In this 
paper, time-scales will be compared to a related "reference 
time^". 



2 We follow the definition of 'Shapiro & Mar chanH fT973) . This refer- 
ence time has a sligh tly different normalization than the relaxation time. A 
IBahcall&WoiafT976 ) cusp forms on a time-scale considerably shorter than 
To, see also equation 



To; 



3n 



16 (Gm)2nolnA 
A/. 



: 22 Gyr 



3 X lO^Mp 



5/4 



(9) 



where the dependence of the Coulomb logarithm on Af, was 
neglected in the last line. 

In the vicinity of a MBH (r < rh), orbits are assumed to 
be Keplerian, and the number of stars N{E)dE in the energy 
range {E, E + dE) is related to the DP by 



N{E) = TT^V2{GM,fE-^^^f{E). 



(10) 



iBahcall & Wolj d 1 9761) showed that the DP for bound stars 
in steady state is (see also lSari & Goldreichl2006l) 



f{E) cx E 



1/4 



N{E) oc E-^/\ 



(11) 



If iVo is the total number of stars that are bound to the MBH 
with energies larger than Vq, then number of stars in the inter- 
val (In E,\nE + d\n E) is given by 



dN 



= -No{ — 



-5/4 



d\nE 4 
The density of stars is related to the DP by 

rGM./r 



(12) 



dEf{E) 

-OO 



GAf. 



E 



(13) 



The critical scale for the semi-major axis of binaries at a 
distance tq of the MBH is set by the hard/soft boundary in 
equation ([3j. Setting ^ = 1 and estimating the squared veloc- 
ity dispersion to be v'\- = GM,/ro, this gives 



ah = 0.07 AU 



M. 



3 X 106 Af© 



-1/2 



(14) 



Larger MBHs live in systems with larger velocity dispersions, 
and hence require tighter binaries in order for them to survive 
in their vicinity. 

We end this section by summarizing some of the assump- 
tions made throughout the paper The assumptions are further 
discussed at several places in the text. 

All stars are single mass, Newtonian point particles, orbit- 
ing a static MBH that is much more massive than the stars. 
Binary stars interact with single stars, the latter having a fixed 
distribution as in equation ( fTTI ). Binary-binary interactions 
are neglected, as is dynamical friction. We consider only the 
evolution of the binaries that are bound to the MBH and as- 
sume a reservoir of unbound binaries that do not evolve. The 
properties of the reservoir are related to those of the MBH 
through the relation between mass and velocity dispersion (01). 
Evolution of energy and angular momentum is entirely due to 
uncorrected, weak two-body interactions, such that recoil af- 
ter an encounter of a binary with a single star is neglected, as 
are secular effects such as resonant relaxation. The internal 
semi-major axis of the binaries evolves only as a result of in- 
teractions with single stars, and we use cross-sections found 
in three-body calculations performed by other authors. 
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3. EVOLUTION OF ENERGY AND ANGULAR MOMENTUM OF THE 

BINARIES 

In this subsection we consider the evolution of the energy 
E and angular momentum J of the center of mass of the bi- 
naries. If N{E, J)dEdJ is the number of stars in a volume of 
size dEdJ around a point (£', J), the Fokker-Planck equation 
can be written as 



dN 
'dt 



d 1 



52 



dEdJ 



[DejN] 



(15) 



The notation can be simplified by introducing dimension- 
less quantities 



t/To; v^E/vl; j = J/ME), 



(16) 



where Jc{E) — GM,{2E) is the angular momentum of 
a circular orbit of energy E. Equation ( fTSl l then becomes 



dN 

s7 



dj ^ 



dydj 



D,N] 



1 92 



(17) 



This partial differential equation is in general non-linear: 
the diffusion coefficients depend on the stellar distribution, 
and this distribution is in turn affected by the diffusion coeffi- 
cients. 

Direct numerical integration of the non-linear Fokker- 
Planck equations is possible, but not straightforward (see 
ICohn & Kulsrudiri978i) . The system considered here is more 
compUcated than the Fokker-Planck equations here describe, 
since evolution of the internal parameters of the binary (semi- 
major axis, exchange probability) is also followed. Instead 
of solving the full problem, we assume simple forms for the 
diffusion coefficients which are determined by a fixed popu- 
lation of single stars that has a iBahcall & Wolfl (11976) cusp 
(see ®. We then follow the evolution of binary stars as they 
interact with the single stars. This leads to a considerable sim- 
plification of the equations, which now become linear. 

Another simplification we will make is that three-body ef- 
fects on {E, J) evolution are neglected. After a strong en- 
counter, a hard binary recoils with a velocity ^ y/GmJa, 
which can lead to its ejection. However, for every such ejec- 
tion, there are of order '--^ In A ejections at similar velocities 
due to the cumulative effect of weak encounters^. The energy 
given in strong encounters to single stars is also smaller by or- 
der 1/ In A than that due to weak interactions, so our assump- 
tion of a static background of is not affected by the negligence 
of strong encounters. This is analogous to the role of strong 
single-single encounters on the dynamics, which was shown 

We only follow the bound binaiies. For each binary ejected with veloc- 
ities much larger than the escape velocity, several binaries will be ejected at 
lower velocities, so the effect of strong encounters on the binary population 
is in fact even smaller than suggested. 



to be unimportant lFreitag et al.l(l2006h . The assumption is fur- 
ther justified by the fact that very little hardening of binaries 
is observed in our simulations (see figure[8]l. 

3.1. Diffusion in energy 

The random motion in energy space due to two-body in- 
teractions leads to energy diffusion. We make the simple 
assumption that the single stars are distributed according to 
the solution found by Bahcall & Wolf (1976), /^(y) oc y^/"^. 
The energy diffusion term D yy was estimated analytically by 
lLightman & Shapirol (1 19771) to be 



85y' 



9/4 



(18) 



(see their equation 51b for p = 1/4). Intuitively, the power 
9/4 can be estimated by noting that the relaxation time is 
inversely proportional to the DF, tr{y) oc y^^^'^, and thus 
Dyy ~ yVUiy) ^ y^y^l^ ^ y^/\ 
A detailed analysis shows tha t in steady sta t e, the 

two-body drift ter m vanishes (iBahcall & Wolll 119761: 

iLightman & Shapirol Il977h . This result was obtained for 
single mass stars. Since our binaries are twice as massive as 
the single stars, there will be some drift in this case because 
dynamical friction becomes more important. Nevertheless, 
we neglect this effect for simplicity, setting Dy = 0. 

3.2. Diffusion in angular momentum 

The angular momentum diff usion term Djj was estimated 
analytically by ILightman & Sh apiro ( 197^ to be 

Dj, = 6.62/1/4 (19) 

(see their equations [55] and [56] for p = 1 /4). Intuitively, the 
relaxation time is the time-scale for the angular momentum to 
change by order of the circular angular momentum, or for j 
to change by order unity. 

The angular momentum drift term due to two-body interac- 
tions is related to the angular momentum diffusion term by 



D 



■J3 



2j 



(20) 



( ILightman & Shapiroll 19771) . 



4. INTERNAL EVOLUTION OF THE BINARIES 

As the orbit of the center of mass of the binary with re- 
spect to the MBH changes due to encounters with other stars, 
the characteristics of the binary itself also change as a re- 
sult of such interactions. In the model used in this paper, the 
most important quantity that evolves due to interactions with 
stars is the semi-major axis of the binaries. For hard bina- 
ries (^ 3> 1, see equation [3]), the binary typically shrinks in 
interactions with single stars. The binary thus becomes even 
harder, and can potentially survive at closer distances with re- 
spect to the MBH. Soft binaries (i^ <C 1) become softer, and 
will eventually be ionized. In addition to the change in the 
semi-major axis, one of the components of a binary may swap 
with the incoming single star during the interaction. The prob- 
ability that an exchange interaction has occurred is followed 
during the evolution. 

4.1. Evolution of the semi-major axes of the binaries 

iHeggie & Hii? (1 19931) find the probability P(A|i;) that a bi- 
nary with initial energy e > that interacts with a star with 
dimensionless velocity 
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V3 

Vc 



(21) 



has a relative change in velocity 



A = , (22) 

e 

where e' is the new energy of the binary (hardening implies 
A > 0, softening implies A < 0; the convention used here is 
that bound binaries have positive energy e). In equation (ISTT l. 
1)3 is the velocity of the incoming star in the center of mass of 
the system, and the critical velocity Vc is the minimal velocity 
required to disrupt the binary, i.e.. 



1 



(23) 



V m 

(see lHeggie & Huj|1993L just before equation 2). 

The dimensionless probability P{A\v) is related to the dif- 
ferential cross-section for relative energy changes A by 



(24) 



7ra2 dA 

This probability function was evaluated throug h numerical 
three-body integrations bv lHeggie & Hull (Il993h . and a fit to 
the results of those simulation is given in equations (10, 35, 
36, 37, 38) of that paper 

The rate at which binaries are scattered from energy e to 
e' depends on the DF of the densities and velocities of the 
stars around it. For the model of a galactic nucleus assumed 
here, the DF is given by f{E) ^ AE^^^ = A[GM,/r - 
u|/2]^/^ (see equation fTTTi. The DF is normalized such that 
/ (Pvf{E) = n{r), where n(r) is the number density of (sin- 
gle) stars. The rate at which the energy of a binary changes 
from e to e' is then 



47rA /•V2GM./r 



GM, 1 o 

T^vi 

r 2 ^ 



da 

""'dA- 
(25) 



The upper limit of the integral is dictated by 0) = 

0; for the lower limit, see below. Using Eq. (l2Tl i to make the 
velocities dimensionless, and eliminating da/dA in favor of 
P{A\v) gives 



(26) 



f.V2GM./(ri.=) 

X / dvv 



GM, 1 , 

V 

rv'i 2 



1/4 



P{A\v). 



This expression can be further simplified by using equation 
( |23] ), and by defining 



(27) 



GM,m 

The expression for Q which is then obtained is 



Q(e, e') = 3^/^Tr^AG^m'/^e-''/^Rix, x'), (28) 



where 

R{x,x') 





■ 1 


1 2 


/ dvv 




V 




3x 


2 



1/4 



P{A\v). (29) 



If the binary hardens, the lower limit of the integral is zero. 
If the binary softens (A < 0), the maximal amount of energy 
that the binary can gain is v"^, because the intruder has not 
more kinetic energy. As a result, the lower limit of the integral 
is given by 



[max(0, -A)] 



1/2 



(30) 



The upper limit of the integral follows from the assumption 
that the single stars are bound to the MBH. 
To find A, we normalize such that at r = r/, = GM,/vg 



ryj2GM./r 

no — AnA / v'^dv 
Jo 



GM, v^ 



1/4 



8.88^1;, 



7/2 
I 

(31) 



so that 



A = ^v-"\ 



(32) 



Combining ( [28T l and ( |32] | then gives 



-3/4 



mvp. 



tR{x, x'). 



(33) 

Rather than considering all the possible energy changes for 
a given value of x, we make a Fokker-Planck approximation 
(see also appendix A) in which only the first and the second 
moments of Q are considered, so 



(A"i?) : 



dAA"i? 



dAA" 



(34) 





■ 1 


1 2" 


/ dvv 








3x 


2". 



1/4 



P{A\v) 



The most negative relative energy change for the binary that is 
considered is A = — 1, because the binary becomes unbound 
for larger changes. Also, the most negative energy must be 
smaller than the energy of the stars moving at the escape ve- 
locity at a given distance from the MBH. This leads to the 
lower limit 



(35) 



3a; 



of the integral in equation (|34| ). Numerical integration and 
fitting then yields the following results for 0.03 < x < 300: 



(00,01,02,03) = (-1.03,-1.23,0.0456,-0.00729) 



(36) 



ln((Ai?)) (Inx)'' ((Ai?) > 0) 

60 = -2.01; 61 = -0.0113; 62 



0.353; 63 = 0.0312. 

(37) 
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ln(-(Ai?» =^c,(lnxr ((Ai?)<0); 

Co = -5 



i=0 



ci = -6.99; C2 = -2.44; cg = -0.356. 

(38) 

The results from the integration and the fits are displayed in 
figure ([T])- 



100 
10 
1 

0.1 
0.01 
0.001 



A 



<AR> 
-<AR> 
<A^R> 
Fit square 
Fit pos 
Fit neg 



le-04 

0.01 0.1 1 10 

x=er/(GMm) 

Fi g. 1 . — The Fokker-Planck coefficients in equation i34\ and the fits 1361 



1 



1 



The rate of exchange interactions per binary is then 



= 53 
or in units of Tq, 

^cxTq = 3 

In these expressions 

/2/3x 



mvp. 



1/4 



10 



1 1 

3x ~ 2' 



1/4 

^0 



1/4 



(41) 



(42) 



(43) 



4 2 81 6 
-v^ H 

a 160 



(44) 

For this derivation similar arguments and definitions were 
used as in ^4.11 A fit to numerical integration of equation 
100 1000 yields 



i?ox = 0.172x-^/4(a;2 + 0.04)-^/^ (45) 
The function (l44l i and the fit ( l45l ) are plotted in figure (|2]i. 



4.2. Exchange interactions 

After an interaction between a binary and a single star that 
does not disrupt the binary, the new binary may be composed 
out of the same stars it started with before the interaction, or 
one of the the original components may be exchanged for the 
single star. This is of importance, since it provides a mecha- 
nism that is likely to exchange main sequence stars for com- 
pact remnants, which may lead to the formation of X-ray bi- 
naries. The relevance of this mechanism is highlighted by 
globular clusters, in which the rate of interactions of hard 
binaries correlates tightly with the number of X-ray sources 
(Poolev et al. 2003; Fregeau 2008) and cataclysmic variables 
dPoolev & H ut 2006). 

The cross-section for an e xchange in teraction to happen for 
a soft binary is found to be (iHutll 19831 see equation 5.3) 

= ^Tr{Gmfmh-^v-^; {v > 1). (39) 
81 

The dimensionless velocity v was defined in equation (Zl[ . 
The exch ange cross-section for a hard binary may be esti- 
mated by dHeggieetalJI 19961) 

a'f^'' ^ a^iGmfm^e-^v-^; = 4.5 (i; < 1). (40) 

In these equations is the relative velocity, and v is defined 
as in equation (l2Tl i. The numerical factor a was found by 
summation of the average of the single mass cross-sect ions for 
direct and resonant exchange in table (1) from Heg gie et alJ 
(11996). 

These cross-sections can be interpolated to find the general 
cross-section for any velocity. 




le-06 ' 

0.001 0.01 0.1 1 10 100 1000 10000 

x=er/(GMm) 

Fig. 2. — The function i?ox in equation (44) as a function of dimensionless 
binary energy x = tr / {GM,m). 

In the simulation, the probability that a given binary has 
experienced an exchange interaction is estimated as follows. 
During every time-step dti at a time ti, the probability that 
there was no exchange interaction is 

i^"°^" = l-rexdi^. (46) 

The probability that a binary has experienced an exchange 
interaction after a time t — Y^idti is then 

= 1 - I{,P^° (47) 
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5. DESTRUCTION OF BINARIES 

There are two major ways in which binaries can be de- 
stroyed. One channel for binary destruction is that (soft) bi- 
naries interact with field stars, and become gradually softer 
If the binding energy is less than some critical value, the 
two components in the binary are no longer considered to be 
bound to each other We will refer to this proc ess as "ioniza- 
tion" following common nomenclature (e.g. iHeggie & Hull 
I2OO3I) . 

The other mechanism that destructs binaries is when a bi- 
nary star comes at peri-apse so close to the MBH that the tidal 
force exerted by the MBH exceeds the internal force of the 
binary components. This happens when the peri-apse of the 
binary is smaller than the tidal radius. 



< ''t = U:r 



2m 



1/3 



(48) 



This happens typically on highly eccentric orbits, so that the 
angular momentum of the binary with respect to the MBH is 
smaller than the angular momentum of the loss-cone. 



J < Ju 



^{l + e)GM,rt. 



(49) 



Theoretically, a binary can also be disrupted when its en- 
ergy with respect to the MBH becomes larger than the energy 
associated with a circular orbit with radius r = rt, but this 
is extremely unlikely, since the time-scale to rea ch J < Jic 
is much shorter (.Frank & Rees 1976: Lightman & Shapirol 
Il977h . 

5.1. The loss cone 

When a binary has J < Jic, this does not necessarily imply 
that it will be disrupted, since the angular momentum can be 
scattered while the star is going to peri-apse. The event rate 
of los s-cone disruptions was first studied by (iFrank & ReesI 
Il976t Lightman & Shapiro 1977 ), and many later papers 
(e.g., i Bahcall & WolJ 119771: iMagorrian & Tremaind [T99l 
ISyer & Ulmer 1999.) . A distinction is made between a "full 
loss-cone" regime, where the orbital angular momentum 
change is much larger than the size of the loss-cone, and an 
"empty loss-cone" regime, where this change is much less 
than the size of the loss-cone. 

The change in angular momentum per orbit is AJp w 
^ PjDjj. If AJp ^ Jic, the binary can jump in and out 
of the loss-cone on a dynamical time, and the loss-cone in 
this regime is always full. The dimensionless rate (in units of 
To) at which a binary with semi-major axis a and energy E is 
disrupted is then approximately 



7fun 



To 
P ' 



(50) 



The demarcation between the empty and full loss-cone 
regime can be found by setting A Jp = Jic- Evaluating this at 
the radius of influence gives (see equations |f~ 



afuii 0.1 AU 



M. 



3 X IQ^Mo 



-5/6 



(51) 



If a > flfuii, then the binary is in the empty loss-cone regime 
throughout the entire cusp; the loss-cone becomes full only 
for stars at radii > r^- In contrast, if a < Ofuu, it is in the 



empty loss-cone regime close to the MBH, and in the full loss- 
cone regime in regions further away from the cusp. 

In the empty loss-cone regime, any binary with J < Ji^ 
is immediately disrupted. The dimensionless rate at which a 
given bina ry enters the loss-cone is of order 1/tr, or, more 
precisely dLightman & Shapirolll977h 



Tn 



Tcmpty " 



tr{E) \n{Jc/Jlc 



E 



1/4 



ln( Jc/J/c) 



(52) 



Since the loss-cone introduces an scale to the system, it 
must be trea ted carefully in the MC simulations. See ^and 
IShapiro & Marchant ( 1978) for details. 

5.2. Ionization 

The rate at which soft binaries lose energy due to the cu- 
mulative effect of many weak encounters can be estimated 
as follows. The change 5v of one of the components of 
the binary in an encounter with a field star at a distance b 
is 5v ~ Gm/bvdisp- On average these encounters cancel, 
(Sv) = 0, but the root mean square of the velocity increases, 
so that ((5e) ^ 6v^. The rate at which single stars pass one 
of the binary components at a distance between (6, b + db) is 
nsVdispbdb, so that the rate at which the energy of the binary 
increases is e ~ {GmYus In A/udisp. The ionization time for 
soft binaries tjon ~ e/e is then 



tu 



'^'disp 



GusTna In A 



(53) 



In the context of ^ the dimensionless ionization rate may 
be estimated by the rate at which binaries change their semi- 
major axis (harden/soften) in three-body interactions at a rate 
(equations !^ [34l i 



(A2i?(l))' 



(54) 



Here (A'^i?(l)) « 0.36, and (A^i?(x)) is approximately pro- 
portional to {E/ef^'^. The rate 7a can be interpreted as the 
ionization rate if x ^ 1, or as the hardening rate when x ^ 1. 

6. DYNAMICAL RATES 

The dynamical processes described in the previous sections 
can be characterized by their rates, which provide insight in 
the results of the simulations presented below. Here the dy- 
namical rates are summarized. All rates are given in units 
of [I/Tq]. Once the rates are known, a simple model can be 
made that describes the DP and the binary disruption rates. 
This model is in qualitative agreement with the results from 
the MC simulations and gives some insight in the dynamical 
processes. 

6.1. Rates 

The rate for an energy change of order unity due to two- 
body interactions is (Eq. [TsT i 



_ DeeTq E 



£2 

only weakly dependent on energy. 



1/4 



(55) 
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While hardening or softening, binaries have a probability 
to have an exchange interaction during every three-body en- 
counter The characteristic rate for these is 



7o 



0.5 



10 



1/4 



i?cx(l)' 



(56) 



In this expression, -Rcx(l) ~ 0.17, and 7ox ~ const for very 
soft binaries, and 7cx ^ x^^'^ for very hard binaries. 

The hardening/softening rates was given in equation (|54] |. 
and the effective loss-cone rate is (see equations iSOl and 1521 ) 



7ic = min[7fuii(i;),7ompty(£')] 



(57) 



Let us consider what is likely to happen to a star that enters 
the cusp on a marginally bound orbit, E k, v^. If the binary 
is very soft (e <C mv^), then the ionization rate 7ion is much 
larger than any of the other rates, and the binary will quickly 
be ionized due to interactions with the other stars. If the bi- 
nary is hard, it is more likely to be tidally disrupted than to be 
ionized. However, the tidal disruption rate is smaller than the 
energy diffusion rate, and it is more probable that the binary 
diffuses to higher energies with respect to the MBH. While 
this is happening, the binary can harden somewhat, but it is 
important to stress that the hardening rate is always smaller 
than the diffusion rate. This implies that when a binary moves 
to higher energies, it becomes usually softer: the semi-major 
axis of the binary may decrease somewhat, but the velocity 
dispersion of the stars in its vicinity grows more rapidly than 
the circular velocity of the binary increases. 

As the binaries diffuse to higher energies, some are tidally 
disrupted. However, since in the empty loss-cone regime the 
tidal disruption rate is always smaller than the energy diffu- 
sion rate, 7cmpty 7diff , this does not lead to a very signif- 
icant depletion of binary stars compared to the situation that 
tidal disruptions are neglected; with respect to this process, 
binaries could exist at all energies"^. This is not true for ion- 
ization of binaries: since hardening occurs at a smaller rate 
than energy diffusion, binaries with any energy a will even- 
tually reach an energy Eg (a) with respect to the MBH where 
they become soft, and for energies E 3> Es{a), 7ion ^ 7diff ■ 
As a result of this there will be an exponential depletion of 
binary stars compared to the situation where ionization is ne- 
glected. 

The rate of exchange interactions for a given binary is com- 
parable to the hardening rate. This implies that 7ox < 7diff , 
and therefore most binaries move first to higher energies and 
are then ionized before they have experienced an exchange in- 
teraction. The probability for an exchange interaction is there- 
fore low. This general result is confirmed in the MC simula- 
tions below. 

6.2. Parametrized model 

The insight obtained from the dynamical rates presented 
here can be used to make a simple analytical model of the 
DF and the disruption rates. We present a model with several 
added parameters to fit the results of the MC model as well as 
possible. The model does not accurately treat the full three di- 
mensional dynamics, but does capture a number of important 
features, and helps to understand the underlying dynamics. 

As argued before, loss-cone effects cannot affect the shape 
of the DF by much for hard binaries, whereas an exponential 

** This is no longer t rue wh en resonant relaxation becomes important, see 
IHopman & Alexandeii )2006al) 



cut-off is expected for soft binaries. The DF can therefore be 
written as 



/(i?)cxi?i/4exp [-E/iw^vD] , 



(58) 



where Wg is a numerical parameter which we chose to fit the 
data best. The normalized number of stars per log E is then 



1 dN 
N^dlnE 



E 



-5/4 



= 7 — exp [-E/iws^v^) 



(59) 



where A'o is the total number of binaries that are bound to the 
MBH (which is much smaller than the total number of single 
stars bound to the MBH). In these expressions, we use the 
initial hardness parameter ^ the binary had when it entered 
the cusp. In reality, the internal energy of the binary evolves. 
However, since 7a ^ 7diff (equationsl54landl55]). the internal 
energy of the binary in fact changes very little as it diffuses 
through the cusp. 
The tidal disruption rate per log _E as a function of energy 



5wt f E 



-5/4 



exp [~E/[w,S^vl)\ 



lie 



lie + la 



(60) 



where wt is a parameter to better fit the results, and a factor 
lie/ [lie + 7a) accounts for the fact that the binary can also be 
disrupted by ionization, decreasing the probability that it will 
be tidally disrupted. 

The ionization rates can be estimated as follows. For bina- 
ries with < E, the rate is 7a. Binaries with > E 
are hard on circular orbits, but for a thermal distribution of ec- 
centricities, a fraction of stars ~ E/ {wg^.v^) spends a fraction 

of time ^ \E/ {wsS^Vq)\ ^^'^ in regions where GM,/r > S^Vq, 
where Wg is the same parameter as in equation ( l59b . We there- 
fore estimate the ionization rate to be 



5Wi 



-5/4 



exp [-E/{ws^vl) 



X mm 



E 



5/2' 



ll 



lie + la 



(61) 



where Wi is the third and last parameter to better fit the results. 

Fitting by eye, we found that a reasonably good approxima- 
tion for the numerical results is given by taking the parameters 
to be {ws,Wt,Wi) — (8, 1,5) (see also figure|7]below). 

7. MONTE CARLO SCHEME 

The evolution in the previous sections been described in 
the terms of Fokker-Planck equations. Such equations can 
be solved using a MC method. The relation between diffu- 
sion equations and MC methods is reviewed in the appendix. 
An advantage of MC methods is that it is relatively straight- 
forward to add effects such as the probability for exchange 
interactions 04.2I ) to such a method. 

In a MC scheme, the DF is represented by a number of bi- 
naries that represents a realization of the actual distribution. 
This number of simulated binaries can be larger or smaller 
than the actual number of binaries. At every time-step, all 
monitored quantities for a given binary change, some deter- 
ministically (if there is a drift), and some randomly (if there 
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is diffusion). In addition, the binary can be disrupted by loss- 
cone effects or three-body effects, there is a probability that it 
experiences an exchange interaction, it can become unbound, 
and so forth. 

In the MC code, the energy and angular momentum of the 
binary are first updated according to 

E{t + St) = Eit)+xi[DEESt]'^^; (62) 



J{t + St) = J{t) + DjSt + X2 [DjjSt] ^1"^ . (63) 

In these equations, xi ^nd X2 are standard normally 
distributed random variables that have a correlation 
De.j / \/ DeeD.t }■ For details see I Shapiro & MarchanB 
(IT978h . 

e(t + St) = e{t) + D,St + V [D.Jtf'^ , (64) 
with ij) an independent standard normal random variable, and 

in A \ / in A \ mwg / 

(65) 

(see equations [33] — |38T l. 

Once the change in internal energy is calculated, we update 
the probability that there has not yet been an exchange inter- 
action for that binary using equation (|46] |. 

A complication in MC simulations is that often the proba- 
bility distribution varies by orders of magnitude over the range 
of interest. This implies that a very large number of particles 
is required in order to resolve the probability distribution at 
values where it is low, and hence that much more particles 
are present at values where it is high. In the situation that is 
considered here, the number of stars per unit energy, dN / dE, 
changes very steeply near a MBH. For example, for a cusp 
with DF / - EP, dN/dE cx E'^/^+p. This impUes that if 
one were to set up a MC simulation naively, a very large num- 
ber of particles needs to be present at energies of order E ^ 1, 
in order for there to be only a few particles at large energies. 

This problem can be resolved by using a "c loning scheme", 
as introduced bv lShapiro & Marchanll (Il978h . We use a sim- 
plified version of their scheme here, which we now briefly 
summarize. 

All binaries start at small energies i?ncw = 0.5wq, with an- 
gular momentum drawn from an isothermal distribution in the 
range J/Jc = (0, 1), and with a semi-major axis drawn from a 
distribution that depends on the model (see next section). For 
each individual particle an adaptive time-step is computed, 
following the prescription presented in Shapiro & Marchant 
(fT9781) . and additional criteria to ensure that the fractional 
changes in the binary properties are smaller than 0.01. 

As the energy of the binary with respect to the MBH 
changes because of two-body and three-body scattering, it oc- 
casionally diffuses to higher energies. If the binary crosses an 
energy E'cione.n — 10", where n — (0,1,2,3), 9 identical 
binaries are made. All these binaries are tagged to be clones. 
The binaries then proceed to evolve independently. When a 
clone crosses the boundary E'cionc.n at which it was produced, 
the binary stops to exist and is discarded in the analysis. A bi- 
nary can also be annihilated by physical processes. The chan- 
nels for annihilation are disruption by the loss-cone; ioniza- 
tion; crossing a large, maximal energy E'max; and interactions 



Table 1 

Model parameters and GW rates 



Name 


[Mq] [AU] 


Q-max 

[AU] 


5o 


Comments 


I 


3 X 10« 







No J and a evolution 


Ila 


3 X 10^ 0.24 


0.24 


0.3 




lib 


3 X 10'=' 0.073 


0.073 


1 




lie 


3 X 10^ 0.0073 


0.0073 


10 




Ild 


3 X 10^ 7.3 X 10-** 


7.3 X 10-"' 


100 




Ilia 


1 X 10^ 13 


13 


0.3 




Illb 


1 X 10^ 4.0 


4.0 


1 




Ille 


1 X 103 0.4 


0.4 


10 




Illd 


1 X 10^ 0.04 


0.04 


100 




IV 


3 X 10^ 0.01 


0.3 


0.24-7.3 


Favored MBH model 


V 


1 X 10^ 0.01 


10 


0.4-400 


Favored IMBH model 



Parameters for the simulations. Models I and II are given for theoretical purposes; 
note that some of the semi-major axes in those models cannot hold for MS stars. 



that lead to an orbit that is unbound to the MBH, where bi- 
naries are only considered to be "bound" by the MBH if their 
energies are larger than Eb — 0.4tig. We consider a binary to 
be ionized when e/rribE < 10^^. In case of these physical 
annihilations, the process that destroys the binary is recorded, 
so that it is possible to distinguish between ionizations and 
loss-cone captures. If the destroyed particle is a clone, it is 
not replaced. If it is one of the original particles, it is replaced 
by a new particle with energy Encw- In this way, the num- 
ber of particles in the simulations that are not clones remains 
constant by construction. 

Using this method, there are approximately equal num- 
bers of simulated binaries at all energies. There is no over- 
sampling of the low energy region, and there are enough par- 
ticles at high energy regions in order to resolve those. At 
high energies, however, particles have less statistical weight: 
any binary with energy 10" < E/v^ < 10"'^^ has a weight 
w = 10"". This weight is then taken into account when cal- 
culating the DF, disruption rates, distribution of semi-major 
axes, etc. Typically, 1000 particles are used at the start of 
the simulation, which procreate to a final number of about 
5000 (depending on the simulation). The distribution reaches 
steady state on a time-scale of the order of the relaxation time. 
These simulations typically completed in a few hours. 

8. RESULTS 



A number of different models are considered to probe how 
the different effects of binary evolution, loss-cone disruptions, 
and the dependence on the initial conditions modify the result. 
The different models are summarized in table ([TJ. 

Model I considers only energy E evolution. The semi- 
major axis and angular momentum of the binaries does not 
evolve, and the re is no loss-cone. This reproduces the 
Bahcall & Wolf ( 1 19761) distribution to good approximation, 
and is used for comparison with other models. 

Models Ila — lid have a il/. = 3 x IO^Mq MBH, and 
include all the dynamical effects discussed in the paper (such 
as J and a evolution, loss-cone). For each of these models, 
the binaries have initially all the same semi-major axis. This 
makes it possible to compare the dynamical evolution of dif- 
ferent binaries. The semi-major axes are chosen such that the 
hardness ratio equals = Gm/2avQ — (0.3, 1, 10, 100) at 
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the radius of influence (see equation |3]l. We note that these 
values were chosen to study the dynamics, and not to be nec- 
essarily realistic. For example, some of these values imply 
that a < Rq ; these cases can be thought of as representing 
compact remnants, or as for theoretical interest. 

The same assumptions are made for models Ilia — Illd, 
with the only difference being that the mass of the MBH was 
now chosen to be Af, — 10'^ A/©, as is perhaps the situation 
in some globular clusters. 

Models IV and V are the models that are considered to be 
the most realistic representation of the nucleus of a Galaxy 
and of globular clusters, within the limits of our model. They 
again include all the dynamical effects. The initial distribution 
of the binary semi-major axes was chosen to have the simple 
form log a, dN/d In a cx const, which is roughly consistent 
with the log-no rmal period distribution of b inaries in the Solar 
neighborhood jDuquennov & Mavodl 19911) . For the Galactic 
nucleus model, the minimum is Omin = 0.01 AU, and the 
maximum Cmax = 0.3 A U, in accordance with the range of 
semi-major axes used by lYu & Tremaind (l2003h . The min- 
imal semi-major axis is the distance for Solar type stars to 
nearly touch each other. The maximal semi-major axis is cho- 
sen such that ^0 ~ 0.3; binaries wider than that will quickly 
dissolve. For the globular cluster model, the minimal semi- 
major axis is the same, but the maximal semi-major axis is 
a^ax = 10 AU, so that ~ 0.5. 

The figures of models I-V all show on the left hand side the 
case of M, = 3 x lO^Af©, and on the right hand side the case 
of M, = lO^M©. Top figures show models II and III, where 
all binaries start with the same semi-major axis, while bottom 
figures show the favored models. 

8.1. The distribution function 

Figure (|3]l shows the DF f{E), integrated over all angular 
momenta for several different examples. If tidal disruptions 
and ionizations are neglected for distances much larger than 
n, a f{E) oc £;^/^ lBahcaU & Wolj(ll976h profile is obtained, 
which provides also a test for the MC code. This profile is 
modified for cases where the loss-cone and three-body inter- 
actions are included. Binaries become ionized when the reach 
high energies, and there is an exponential cut-off to the DF. 
For the cases where at the radius of influence ^ = 0.3 (mod- 
els Ila and Ilia), the DF quickly drops to zero. It is clear that 
such binaries can in fact not exist in the vicinity of MBHs. For 
smaller semi-major axes (^ > 1, models b-d), binaries can to 
some extend survive on orbits bound to the MBH, and the 
smaller the semi-major axis, the higher the energy to which 
it survives. This is the combined result of the fact that tight 
binaries are less likely to enter the loss-cone of the MBH, and 
they are less easily ionized. The qualitative comparison be- 
tween the MC results and the parametrized analytical model, 
shown for A/, = 3 x 10^ Af©, is good. Comparing the fig- 
ures for M, = 3 X lO^Af© and A/. = 10^ A/©, it can be 
seen that the results are very similar. This is to be expected: 
as discussed in ^ the only process that can strongly modify 
the DF is the annihilation of binaries by ionization. The equa- 
tions determining the rates for this to happen only contain the 
dimensionless quantity ^ which relates the internal energy of 
the binary to the velocity dispersion around the MBH. In the 
figures we scale the energies with the velocity dispersion at 
the radius of influence, so that any dependence on the MBH 
mass is scaled out. The only MBH mass dependent process 
is then disruption by tidal effects: the critical distance where 
the loss-cone becomes full is M, dependent. However, in ab- 



sence of resonant relaxation, as assumed here for simplicity, 
loss-cone effects cannot strongly modify the DF. The binary 
disruption rates will be further discussed in figure (|6]l. 

The bottom figures show the DF for the Galactic center 
(left) and IMBH (right) model. Binaries can exist at much 
higher energies for our "realistic" IMBH model compared to 
the "realistic" GC model. The reason for the difference is 
that for real systems, the stellar radius is an additional length 
scale of importance that starts to play a role. Binaries cannot 
be tighter than their own radius, and in this model we there- 
fore allowed for much harder binaries in the IMBH model 
(C < 400 for IMBH compared to ^ < 7.3 for GC, see table 
[ij. If instead one would consider the evolution of binaries 
of compact remnants, their radii, or the radius where gravi- 
tational wave emission becomes important, are so small that 
there would not be a difference betwe en the IMBH and GC 
models (see also lO'Learv et al]|2008l for the role of stellar 
black holes in galactic nuclei). 

Another way of viewing the effect of ionization on the bi- 
nary DF is shown in figure (|4]i. In this figure, the two dimen- 
sional distribution E^^'^d'^N{E, a)/d\nEd\n a is displayed; 
the factor E^^^ was added to factor out some of the energy 
dependence of the DF (see equation |59]|. Again it is clear that 
at small energies binaries of all types can exist, whereas at 
high energies there is an a-dependent maximum energy that 
is approximately given by Es{a) — WsGm/{2a) — AGm/a. 
Here = 8 is the same parameter is used in the model in 



8.2. The binary fraction 

The number density of stars is related to the DF by equation 
( fTsl l. Let the binary fraction be defined as the fraction 
of the number of binaries in a particular scenario, compared 
to the number of binaries in the absence of a loss-cone and 
three-body interactions. 



_ ribW ^ / <Pvfb{E,J) 
- np{r) Jd^vfpiE) 



(66) 



(the subscript p stands for "point particles"). The actual bi- 
nary fraction is much smaller than this, and equals the fraction 
of binaries at the radius of influence o, times fb{r). Since 
we did not attempt to model the evolution of binaries outside 
the radius of influence, we keep /{, o as a free parameter (so 
fbir = Th) = 1). 

Even if no binaries at all are bound to the MBH beyond a 
certain energy, there is always a finite density at any distance 
from the MBH because of loosely bound, highly eccentric or- 
bits. The density of such orbits in a Kepler potential increases 
as n cx r^^/'^ towards the MBH. Since the density of single 
stars is Ug oc r^^/'', the binary fraction due to highly eccentric 
orbits is 



^fly-by(^)^^./.^ (67) 

Figure Q displays the binary fraction fbir) as a function of 
radius, for models II — V. As in figure (O, one can see that the 
harder the binary, the higher the binary fraction close to the 
MBH. Again the results for binaries of a given hardness are 
very similar when the IMBH and GC cases are compared, for 
the same reason that led to the similarity between the DFs. For 
the realistic GC model, the binary fraction can be substantial 
outside ^ 10% of the radius of influence. Closer to the MBH, 
the profile drops as fb{r) cx r^^"*, implying that only "fly-by 
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Fig. 3. — DF of several types of binaries near MBHs and IMBHs, as a function of energy. Figure (a) shows models II, figure (fc) sliows models III, figure 
(c) sliows model IV, and figure (c) sliows model V. In all figures we also sliow model I in which there are no tidal disruption of binary evolution effects for 
comparison. The DF is strongly affected by the environment near MBHs, which is hostile to binaries. Soft binaries are more easily disrupted than hard binaries; 
the latter can survive to somewhat larger energ ies. In figures (a) the analytical approximation given in equation (58) is also plotted. The straight line is oc i?^/*, 
which is the theoretical lBahcall & Wolj 11973) result when there are no loss-cone and ionizations. 



binaries" contribute to the density profile there. In contrast, 
binaries can exist in substantial amounts in all regions near an 
IMBH, the reason being (again) that much harder binaries are 
allowed for that case. 

8.3. Binary disruption rates 

Figure @ shows the rate A^^"^ (£■ N / dXnEdt) at which 
binaries are destroyed by ionization or tidal disruptions for 
the different models. For soft binaries = 0.3), the ion- 
ization rate quickly becomes higher than the loss-cone rate 
as energy increases. For binaries that are hard when they are 
marginally bound to the MBH, loss-cone effects dominate at 
low energies. These binaries can only ionize on very eccen- 
tric orbits where they spend a small fraction of time in very 
hot regions near the MBH. At higher energies, ionizations al- 
ways dominate over loss-cone disruptions. Since the time- 
scale for ionizations is much smaller than the relaxation time, 
energy diffusion is not sufficiently rapid to replace binaries 
that have been destroyed, and the stars are depleted at high 
E. Ionizations for binaries of a given initial semi-major axis 



is thus localized around a small range in energies, while tidal 
disruptions decrease gradually as the energy grows, and then 
suddenly drops off because ionization effects annihilate all bi- 
naries rapidly. The analytical toy model captures these aspects 
of the disruption rates. 

The situations for M, = 10^ Af© and M, = 3 x lO^Afg 
are quite similar The only place where the difference in mass 
comes into play is that for binaries of a given hardness ^, the 
point where the full loss cone becomes empty is located at a 
higher energy for MBHs of lower mass. The turn-over from 
the full to empty regime is most clear for the case where Af, = 
lO^Af©, ^ = 100 (model Illd). 

In contrast, the difference between the realistic IMBH and 
GC models is again quite clear As argued before, very hard 
binaries cannot exist near the GC MBH, because of the fi- 
nite size of the stellar radius. Loss-cone rates do dominate 
for very low energies, but ionization effects take over quickly 
as the energy increases, and at an energy of E ^ IOOwq, most 
binaries have been ionized. For IMBHs, the majority of the bi- 
naries are so hard that loss-cone effects dominate everywhere. 
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Fig. 4. — Two dimensional distribution E^^'^d? N (E , a)/d\n Ed In a of binaries with arbitrary scaling for the GC model (figure a) and the IMBH model 
(figure b). The smaller the semi-major axis of the binary, the more resilient it is against ionization, and hence the larger the energy with respect to the MBH for 
which it can still exist. The line indicates the boundary E = WsGm/(2a) = AGm/a where binaries approximately become soft. Since binaries are harder for 
the IMBH model, they can exist at much larger radii. Note that the horizontal axis of the figures have different scales. 



In the GC model presented here, the tidal disruption rate of 
binaries in the Galactic center is 7 x 10~^(/b,o/0.05) yr~^, 
where o is the binary fraction of stars far away from the 
MBH. This rate is lower by an order of magnitude than the 
rate found in lYu & Tr emaine (2003). The cause of the differ- 
ence is that we only considered binaries that are bound to the 
MBH. Figure (|6^) shows that the tidal disruption rate contin- 
ues to increase as the energy decreases, since the loss-cone 
in this regime is still empty for binaries. The tidal disruption 
rate of binaries can e ven be much larger than that given in 
lYu & Tr emaind ( l2003i) if there are massive perturbers at a few 
pc from the MB H, or if the potenti al is significantly triaxial at 
those distances (iPerets et al.ll2007h . 

For the IMBH model, we find a tidal disruption rate of 
10"^(/b,o/0.1) yr"^ In contrast to the GC model, this rate 
will not increase much if binaries unbound to the MBH 
are considered, since most binaries are in the full loss-cone 
regime. The ra te we find h ere is much higher than the 
rate estimated in iPfahll (l2005h . who found an event rate of 
^ 10^^(/b,o/0.1) yr~^. The discrepancy can be traced to the 
fact that the models for the ambient cluster were quite differ- 
ent, with our model having a much smal ler re l axatio n time. If 
we were to use the relaxation time of iPfahll (l2005h . we find 
an event rate of 2 x 10~^(/b o/O.l) yr~^, which is in good 
agreement with his results. 

It is evident that the disruption rates near IMBHs are 
highly uncertain, given that there is not much known about 
such systems. It is interesting to note, however, that cur- 
rent models of tidal disruptions, both of binaries and of sin- 
gle stars, typically lead to event rates that are so high that 
the IMBH dis rupts more than its own mass w ithin a Hub- 
ble time (e.g. IWang & Me^ l2004t iHopmanI 120091) . It is 
therefore unlikely that IMBH systems can exist in this form 
for a Hubble time: either the IMBH mass will grow, or the 
system around the MBH will expand and become less dense 
(iMarchant & Shapiro 1980), leading t o lower event rates ; this 
effect is seen in TV-body simulations (iTrenti et al.ll2007h . An 
obvious third possibility is that IMBHs do not exist. 



In figure (|7]i we show the same curves as in figure (|6};), and 
add the analytical approximations from ^6.21 The qualitative 
agreement is good, and the analytical approximation captures 
the main dynamical effects, such as the localized peaks in ion- 
ization, the energy dependence and magnitude of the tidal dis- 
ruption rates, and the cut-off of all rates at high energies. We 
therefore believe that the analytical arguments presented in 
36.21 are essen t ially c orrect. 

Tre nti et aP (l2007h report that tidal disruption rates of bina- 
ries were much lower in their A^-body simul ations than ana- 
lytical predictions suggest. The binaries in T renti et al.l (|2007l) 
were hard compared to the average kinetic energy in the clus- 
ter, but presumably many binaries were soft compared to the 
environment near the radius of influence, where the velocity 
dispersion is higher than average. Figure ^p) shows that for 
soft or marginally hard binaries, the ionization rate is larger 
by an order of magnitud e than the tidal disruption rate, sim- 
ilar to what Trenti et"al] (|2007|) found. This implies that the 
loss-cone rate is smaller by a factor 7ic/(7ic + 7a) < 0.1 
than one would expect if ionizations are not accounted for 
(see equation ll60l ). This factor brings th e analytical model in 
qualitative agreement with the results by lTrenti et al.l (l2007h . 

8.4. Internal binary evolution: semi-major axis and 
exchange rate 

In figure (O, the evolution of the distribution dN / d In a of 
semi-major axes are shown for times t = and t ^ Tq for the 
GC and IMBH models. There has been some hardening of 
binaries, and some of the soft binaries are destroyed, leading 
to a moderate change in the distribution. The relative change 
in the semi-major axis is smaller for harder binaries, which 
have a smaller cross-section. However, as discussed in ^ 
the semi-major axis distribution cannot evolve by very much 
(assuming that there were no very soft binaries at E/vq = 1), 
because the diffusion and destruction rates are larger than the 
rate for hardening. It is not very clear from the figure that 
"Heggie's law", stating, roughly, that hard binaries become 
harder while soft binaries become softer, is of relevance here. 
Hard binaries do not have time to become harder, and when 
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Fig. 5. — Binary fraction as a function of radius. Figure (a) sliows models II, figure (b) shows models III, figure (c) sliows model IV, and figure (d) shows 
model V. As in figure (5), it is clear that hard binaries can survive closer to the M BH. The straight line is oc r''/^, which is the radial dependence of the binary 
fraction if only highly eccentric orbits contribute to the density profile (see equation|67). Note that these binary fractions are scaled such that /i,(r = r^) = 1; if 
the binary fraction of the type considered is lower than 1 , as will typically be the case, the actual binary fraction needs to be multiplied by the required factor In 
particular, for the Galactic center the actual binary fraction will be smaller by a factor < 0.1. In the IMBH model, binaries are present even at very small radii, 
and the r^/^ limit is not reached. 



they move to higher energies, where the densities are higher, 
they tend to become softer so that they are disrupted. Soft 
binaries do become softer, but this happens so rapidly that 
they are quickly ionized, so that they do not count in the final 
statistics. 

Figure (|9|l shows the distribution of probabilities that a bi- 
nary has undergone an exchange interaction for models (II — 
V). Binaries with ^ = 1 at E/vq have the highest probabil- 
ity: harder binaries have smaller cross-sections, while softer 
binaries are more likely to be ionized and do not survive for 
a relaxation time (this is in particular true for even softer bi- 
naries with ^ < 0.3, not shown in this figure). Nevertheless, 
the probability for an exchange interaction of a given binary 
is never much larger then ^ 0.1, and usually much smaller, 
with an exchange probability of pcx ~ 0.03 being typical. 
The probability is in particular much smaller than what one 
would expect from an "nwcr" argument (or from the pre-factor 
in equation ||56l ). The reason for the discrepancy is that bina- 
ries typically do not survive for a relaxation time in the Galac- 



tic center, but are prematurely disrupted by ionizations (even 
hard binaries can move rapidly to higher energies, where they 
become soft), or in the loss-cone. 

9. SUMMARY AND DISCUSSION 

In this paper we analyzed the evolution of single mass bi- 
nary stars in a fixed environment of single stars and a MBH. 
We followed the evolution of the semi-major axis of the bi- 
naries, the energy and angular momentum with respect to the 
MBH, and the probability that a binary has experienced and 
exchange interaction. Channels for disruption of binaries are 
the ionization by a three body interaction, or tidal disruption 
by the MBH. The main goal was to provide a scheme to study 
the resulting dynamics, which we analyzed using a MC code. 

We have compared simulations with M, = 3 x IG^A/q 
and M, = 10'' Mq. Many of the properties of these systems 
are very similar. However, the mass of the MBH discrimi- 
nates the systems in two ways. First, the tidal radius depends 
on M,, and hence the loss-cone for binaries of a given hard- 
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Fig. 6. — Disruption rate per unit In E, per binary that is present in the system. A distinction is made between loss-cone disruptions and ionizations. Figures 
{a, c, e) are for a M, = 3 X lO'' MBH, figures \b, d, f) are for an IMBH. Figures (e) and (/) are the "realistic" models, the others are cases where all the 
binaries start at a given initial semi-major axis, as indicated in the legends. Note that the top figures have a different scale. For all but the softest binaiies, tidal 
disruptions dominate over ionizations at low energies. As the binaries diffuse to higher energies, they become softer with respect to their environment, until 
they reach the point where their internal energy becoines smaller than the kinetic energy of the stars around thein, at which point they are rapidly ionized, and 
ionization rates exceed tidal disruption rates at sufficiently high energies. This process is so rapid that it depletes all binaries and at higher energies, binaries can 
no longer exist; the ionization rate is therefore localized in energy. The realistic models give tidal disruption rates of 7 X 10~^(/i, o/0.05) yr~^ for the GC and 
10~^(/i, o/O-l) yr~^ for IMBHs, where is the binary fraction of stars far away from the MBH. See text for discussion of these rates. 



ness ^ is different for the two systems. This is quantitatively 



expressed in equation dSTT l. Second, assuming the relation be- 
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ruption rate is thus much higher than we find in our analysis. 
It would be of interest to extend the treatment to consider un- 
bo und orbits; thi s can be done within the same MC scheme, 
see lMarchant & Shapiro ( 1979, 1980). As we noted in ^ it is 
in particular important to include global dynamics for IMBH 
systems, which have short relaxation times and cannot exist 
for a Hubble time in the state proposed here. 

A final limitation we stress here is tha t we did 
not include resonant relaxation (e.g., Rauch & Tremai"n3 



1996; Rauch & Ingalls 1998; Hopman & Alexander' ^006at 
Gurkan & Hopman 2007; Eilon et al. 2008; Perets et alj 
20091) . Throughout this paper, we have argued that loss-cone 
effects cannot dramatically change the DF, since energy diffu- 
sion operates on a compar able (somewhat short er) time-scale. 
It was shown by Hopman & Alexander (2006a) that resonant 
relaxation can change this. This relaxation mechanism is ef- 
fective close to MBHs where the potential is close to a Kepler 
potential, so that orbits precess very slowly, leading to strong 
torques between the orbits. Since the resonant relaxation time 
can be much smaller than the non-resonant relaxation time, 
close to MBHs, energy diffusion (which still acts on the non- 
resonant relaxation time) is too slow to replenish binaries that 
are tidally disrupted. We expect that resonant relaxation time 
will modify the results presented here in the region very close 
to the MBH (- 0.01 pc for M, = 3 x IO^Mq). While 
the MC method can be easily adapted to include this effect, 
we did not incorporate it since the details of resonant relax- 
ation are not yet fully understood, in spite of ongoing efforts. 
In particular, it is not well known how the mechanism de- 
pends on back-reaction effects (known as resonant friction, 
see Rauch & Tremaine 1996). 

As outlined in the introduction, there is considerable astro- 
physical motivation to study binaries in the vicinity of MBHs. 
We now consider a few of the ramifications that could be of 
observational relevance here. 

X-ray binaries in the Galactic center — There is an over- 
abundance o f X-ray binaries iri the Galactic center^ compared 
to the bulge dMuno et"ani2005h . Since there is very little evo- 
lution of the internal semi-major axis of the binary (see ^6.11 
and 38.4l i. hardening cannot lead to additional X-ray binaries. 
Similarly, since the binary fraction drops within the radius of 
influence (see ^8.21 ). the migration of X-ray binaries from the 
bulge to the central parsec can also not account for the rela- 
ti vely high number o f X-ray binaries per unit stellar mass. 

iMuno et al] (l2005h suggested that the over-abundance of X- 
ray sources is the result of exchange interactions. These can 
turn neutron stars into binary members, which can then con- 
tinue to produce X-ray binaries. We now compare their anal- 
ysis to the one made here. The mean probability for an ex- 
change interaction of a binary in the GC is (pox) = 0.03. 
Assuming that a fraction /ns = 0.01 of all stars in the GC 
are NSs, we find that the number of NSs that were born as 
single objects, and are members of binaries at the end of the 
simulation can be estimated as 
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Fig. 7.— Disraption rates for the A/. = 3 X lO", § = (10, 100) models 
(thick lines), compared with the analytical approximation (thin lines with 
same colors for given model). 



tween MBH mass and velocity dispersion of the host galaxy 
(equationlH implies that binaries of given hardness ^ need to 
be tighter as increases. Since stars have finite radii, this 
implies that the lower the mass of the MBH, the harder the 
binaries can be, see equation ( fT4l i. 

It is perhaps worthwhile to stress several of the limitations 
of the approach we present in this paper, without attempting 
to provide an exhaustive list. 

The choice of assuming single mass stars obviously limits 
the use of our work, and it is of interest to consider a spec- 
trum of masses. This is considerably more complicated, not 
only since the outcomes of the three-body interactions will 
then have much more parameters, but also since the dynam- 
ics with respect to the MBH will be mass dependent. In the 
current analysis we have also neglected dynamical friction by 
making the assumption that Dy — 03.1l l. Integrating the 
multi-mass Fokker-Planck equations given in lBahcall & Woln 
([1977), we find that for a two-mass population with the mas- 
sive stars being twice as heavy and ten times as rare, the dis- 
tribution of the light stars is not affected and has an cx r^^-^^ 
profile, while the distribution of the heavy stars (represent- 
ing binaries) is slightly steeper with n (x r^^-^. This con- 
firms that the distribution of the single stars is not significantly 
modified by the presence of the binaries, while the binaries 
themselves are driven towards the MBH by dynamical fric- 
tion, although the effect is quite mild. On a qualitative level, 
this implies that the life time of binaries within the stellar cusp 
is further restricted, and they would be ionized more quickly 
(see also the discussion of the exchange probability below). 

Similarly, the inclusion of binary-binary interactions would 
make the analysis much harder. For both these extensions 
of our analysis, we believe that the study presented here will 
capture much of the dynamics and will therefore be useful to 
compare with. This is especially true for binary-binary inter- 
actions, which are very rare due to the small (^^ 1 — 10%) 
binary fractions. 

Another limitation is that we only evolved the stars bound to 
the MBH. Since the loss-cone is empty for all but the tightest 
binaries for the GC model, most of the tidal disruptions come 
in fact from orbits that are not bound to the MBH (see, e.g., 
lYu & Tremainel2003l;|Perets et alj2007l) . The actual tidal dis- 



^^cx,NS ~ fmfb.o / d^r{pc^)fb{r)ns{r) « 20 



/ns /b,0 

0.010.05' 
(68) 

The number of binaries containing an exchanged com- 
pact remnant is lower by a factor 10-100 than estimated in 

^ we restrict our discussion to the inner pc of the Galactic center. At larger 
distances not only does our method not apply, but the relaxation time exceeds 
the Hubble time, implying that there will be very little dynamical evolution. 
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Fig. 8. — The initial and final distribution per In a of the binaries for the GC model (a) and the IMBH model (b). 
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iMuno et al] ( l2005h . and about 25% of all binaries should then burst every year in order to explain the observed transient X- 
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ray sources. The number of binaries containing a NS can also 
be higher if the fraction of NS /ns is higher than assumed 
here, or if the fraction of unbound hard binaries is larger. A 
possible third cause of the relatively low event rate estimated 
here is our simplifying assumption that all stars have identical 
mass. Since compact remnants are likely to be more massive 
than typical main sequence stars, they are possibly more likely 
to be exchanged with one of the binary components than we 
assumed here. At the same time, more massive binaries will 
drift towards the MBH more rapidly, increasing the probabil- 
ity for ionization. It is not a priori obvious that accounting for 
differences in mass is actually conducive to the exchange rate 
of compact remnants. A treatment of binaries with a spectrum 
of masses is outside the scope of this paper, but could be stud- 
ied with a method similar to the one presented here in future 
research. 

Hyper-velocity stars from the Large Magellanic Cloud — 

For a MBH of mass similar to that in the Galactic center, 
the rate of tidal disruptions is dominated by binaries that are 
unbound to the MBH. This is not the case for hard binaries 
around IMBHs; as can be seen from figure dlJl), the loss-cone 
becomes full at energies E/v^ > 1. We can therefore use our 
simulations to estimate the tidal disruption rate near IMBHs. 

Tidal disruptions of tight binaries lead to the produc tion of 
hyper- velocity stars (iHillslll988tlYu & Tremain3l2003l) . Most 
of the observed stars are consistent with being ejected from 
the Galactic center, but the star HE043 7-5439 appears to b e 
too young to have originated there (E delmann et al.ll2005h . 
This has led to speculations that this star has instead been 
ejected from an IMBH in the Large Magellanic Cloud (LMC; 
see, e.g., Edelmann et al. 2005; Gualandris et al. 2005), an 
hypothesis that appears to be supported b y its metallicity 
dBonanos et al.ll2008t IPrzvbilla et all 120081) . An alte rnative 
idea is that HE0437-5439 is a blue straggler (e.g. iPeretsI 
l2008h. 

iGualandris et alj (l2005h performed three-body simulations 



to find the fraction of tidal disruptions by an IMBH that leads 
to the production of a hyper-velocity star with the properties 
of HE0437-5439. They found that for a « 0.1 AU, a frac- 
tion f^^^ = 0.1/^^^ of the tidal disruptions lead to the re- 
quired velocities. From figure (|6}l), we estimate that the dis- 
ruption rate per binary for such ratios is Nq^cPN/ d In Edt ^ 
0.1 Myr^^. This result depends on the boundary condi- 
tions that are assumed, see 38.31 The boundary condi- 
tions here assume that IMBHs follow the relation between 
M, and vq (see equation |4|, consistent with evidence from 
- /V-body simulations for runaw ay mergers in young clusters 
jPortegies Zwart & McMiiian|[2 002). Let f = 0.1 ^ be the 
fraction of binaries with a < 0.1 AU, and = O.OlfJ^^^^ 

the fraction of stars with mass > 8Mq, then the total num- 
ber of relevant binaries is Nq = 103/f'/™>8 = flJ^>^. 
We thus find that the rate at which an IMBH in the LMC 
could eject hypervelocity stars of the type of HE0437-5439 is 

d^iN/No)/dlnEdtf^^^No - 0.01 Myr-i/o'i/ao^'/oT'- 
Assuming an age of 10 Myr for the age of the star, and ac- 
counting for a geometric factor such that 25% of the ejected 
stars move towards the Galaxy, we find that the probability 
that we could observe a star like HE0437-5439 that originated 
in the LMC can be estimated as plmc 0.03/n ^ fnm.f^Y^ 



This e vent rate is higher t han that implied in Gualand ris et alj 
(l2005h : see iPeretsI (120081) for a discussion of how their ejec- 
tion rates can be converted into detection rates. We conclude 
that the hypothesis of an origin in the LMC cannot be strongly 
rejected, although at the same time such an origin appears not 
to be very likely. 

I thank Atakan Giirkan, Douglas Heggie, Yuri Levin, and 
Hagai Perets and for very stimulating discussions and advise, 
and the referee for a very helpful report. This work was sup- 
ported by a Veni fellowship from the Netherlands Organiza- 
tion for Scientific Research (NWO). 



APPENDIX 

THE FOKKER-PLANCK APPROXIMATION 

At several places in this paper, we have used Fokker- Planck approximation s, and MC methods to solve these equations. In this 
appendix we summarize several (well known, see e.g. IChandrasekharifT944l) aspects of the Fokker-Planck approximation, and 
derive that these equations indeed describe the distribution of MC simulations. 

Consider a general quantity x of a star, which changes due to dynamical interactions (this x has no relation to earlier use of 
X in the paper). If Ax) is the probability for a change of size Aa; of a star at x, then the rate at which the number of stars 
n{x)dx in the range {x, x + dx) changes is given by the master equation 



dn{x, t) 
dt 



J dAxnix - A., tM^- Ax. Ax) - nix, t) J dAx^^, Ax 



(Al) 



The first term of the master equation gives the rate at which stars flow into the bin {x,x + dx), and the second term the rate at 
which they flow out of this bin. The probability distribution ^'(a::, Ax) is determined by the "micro physics", such as two body 
interactions, or interactions of a binary with a single star. 

The master equation is an integro-differential equation. In many situations in physics, the behavior of a system is determined 
mainly by the sum of many small changes in a quantity, rather than rare but strong encounters. In such situations, it is a good 
approximation to make a second order Taylor approximation of the master equation ( lAll ). Approximating 



•^{x - Ax, Ax)n{x - Ax, t) 
and defining 



'^{x,Ax)n{x,t) 



d 

Ax— I'i'ix, Ax)ti(x, t)] 
ox 



1 92 

-Ax^-^[-^{x,Ax)n{x,t)] 



(A2) 



D, 



dAx<i/{x,Ax)A^ 



the Fokker-Planck equation can be written as 



(A3) 
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dn{x, t) 

at 



[Di{x)n{x,t)] + [D2{x)n{x,t)] 



(A4) 



From equation ( IA3I ) it is clear that the diffusion coefficients are just the moments of the scatter probability 'i. 

The first term in the Fokker-Planck equation clearly implies a drift of stars in x. In a MC simulation, for a time-step dt there is 
a deterministic step of size Di{x)dt. 

The second term represents a step in x of size y/ D2{x)dt, but with a random sign (with equal probability to be positive and 
negative). To show this, we derive the diffusion equation for an un-biased one dimensional random walk with varying step size. 

Let the step-size at each time interval Ai be Ax, where Aa; = /S.x{x) is a function of the quantity x of the star First calculate 
the flux T{x^ t) of stars through x at time t, i.e., the net rate at which stars flow from positions smaller than x to positions larger 
than x. To second order in Aa;, the smallest position from which a star can still flow through x is then 



A_ = Aa;(x - Aa;) 



d/S.x . , , . dAa; 
Aa;(a;) ; — Aa;(xJ = Aa; ; — Aa;. 

dx 



dx 



Similarly, the largest position from which particles can step over x is 

d/\x 



A+ = Aa;(a; + Aa;) w Aa;(x) + — - — Aa;(a;) = Aa; + 



dAa; 



-Aa;. 



(A5) 



(A6) 



dx dx 

For an unbiased random walk, half of the particles step right, and half of the particles step left. If at time t the number of 



particles in the interval (a;, x + dx) is given by n{x, t)dx, the flux is then 

1 r 1 
JF(a;,i) = — — / dxn(x,t) — —— 
^ ^ 2At ^ ^ 2At 



dxn{x, t). 



(A7) 



Using J^^ ^ dxf{x) w hf{x + h/2) w f[x)h + ■^^,the flux equals up to second order to 



1 A 1 A 

^ ' ' 2At ^ 2 ' 2At ^ 2 ' ^ 



Defining the diffusion coefficient 



and using the continuity equation 



this leads to the diffusion equation 



1 

2At 
1 

~2At 
_ld_ 
~'2dx 



n{x, t) - 

n{x, t) 
Aa;2 



1 dn 
2'dx 

1 dn 
~ 2'dx 



Ax 
Ax 



Ax - 
Ax 



dAx 
dx 
dAx 



dx 



Ax 
Ax 



At 



n{x, t) 



D2{x) 



dn{x, t) 
dt 



dn{x,t) 1 92 



[Ax(a;)]^ 



At 



dx 



T{x,t), 



[D2{x)n{x,t)] 



(A8) 



(A9) 



(AlO) 



(All) 



dt 2 dx^ 

Since this is the second term in the Fokker-Planck equation (IA4I) . it follows that the Fokker-Planck equation can be solved by 
a MC method in which stars make at each time step a deterministic step of size Di {x)dt and a step of magnitude ^ D2{x)dt and 
random sign. 
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